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Abstract 

An extension of non-deterministic processes driven by the random telegraph signal 
is introduced in the framework of piecewise deterministic Markov processes [9] , in- 
cluding a broader category of random systems. The corresponding Liouville-Master 
Equation is established and the upwind method is applied to numerical calculation 
of the distribution function. The convergence of the numerical solution is proved un- 
der an appropriate Courant-Friedrichs-Lewy condition. The same condition preserve 
the non-decreasing property of the calculated distribution function. Some numerical 
tests are presented. 
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1 Introduction 



Modeling non-deterministic physical systems is a modern subject of research 
that covers many different disciplines. For models that use an a priori source 
of randomness, the Gaussian white noise is the most common ingredient. Nev- 
ertheless there is a growing interest in searching alternatives to the white noise. 
There are two main characteristics, often considered unrealistic for modeling, 
that motivate these researches: the white noise induces (i) infinitely many ran- 
domness for each small interval of time and (ii) unlimited fast fluctuations, 
although with an exponentially decreasing probability. The random telegraph 
signal (named also binary or dichotomous process) is the most elementary 
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assumption that does not suffer of tliese patliologies: it is a signal switching 
randomly between two values, with a determined law for switching times. As 
examples of application to science and engineering, we mention: anomalous dif- 
fusion [6], reaction-diffusion [12], scattering of radiation [13], biological disper- 
sal [11], for diffusive processes, and non-MaxweUian equihbrium [5,8,10,15,19], 
diagnostic technique for semiconductor lasers [14], filtered telegraph signal ^ 
[18,14], harmonic oscillators [16], for systems having an equilibrium. 

For an introduction, let us consider a dissipative process subject to a noised 
input: 



where X{t) is a dynamical variable and ^{t) represents the noise. When ^{t) 
is the Gaussian white noise, X{t) is an Ornstein-Uhlenbeck process, and the 
associated equation dtp{x,t) — dx{xp{x,t)) -\- (l/2)9^p(x, i), describing the 
density distribution function p{x,t), results by the well known Fokker-Planck 
equation [2]. When ^{t) is taken as the random telegraph signal, the same 
equation describes a filtered random telegraph process [18,13] or, equivalently, 
a Langevin equation [5,8] subject to a dichotomous noise, ^{t) alternately 
takes on values ±1, with an exponential (or Poisson) switch probability den- 
sity function of the form: ^e~'^'^, where fi~^ is the average time (r) between 
switches. With this assumption, we see that the Eq. (1) can be integrated as 
an ordinary differential equation, provided that not any switching event hap- 
pen inside the integration interval. Therefore, with the exceptions of switch- 
ing times, the process is deterministic, composed of pieces of increasing and 
decreasing exponentials. Anyway, the whole resulting process X{t) is not de- 
terministic, it represents a random sample path in a probability space, and 
the statistical properties can be described using the characteristic functional 
method [17,8,7,14] or the associated probability densities distribution p^{x, t), 
governed by a Liouville-Master Equation ^ of the following form [20,14,5]: 



This kind of process is a case of a more general class of piecewise-deterministic 
processes that were formalized by Davis [1,9]. Following this framework, in 

^ The author acknowledge K.D. Jakeman and E. Ridley for having advised of 
Pawula's works. 

^ In general, equations for density probability of random processes are derived from 
the Chapman-Kolmogorov equation. As discussed in Ref. [2] the same equation be- 
comes a Liouville equation, in absence of randomness, and a Master Equation, if only 
jump processes are involved. We use both terms in order to stress the deterministic 
and the random character of the processes considered here. 
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Sect. 2 we define an extension of Eq. (1), generating continuous piecewise de- 
terministic processes, and establish the corresponding Liouville-Master Equa- 
tion, fike Eq. (2). Sect. 3 is devoted to investigate the properties of the upwind 
method used to solve the Liouville-Master Equation. In Sect. 4 the numeri- 
cal solution of a dissipative model is plotted and compared to Monte Carlo's 
histograms of the same process. 



2 Continuous piecewise-deterministic Markov processes 

The process described by Eq. (1), rather than a single equation driven by 
a binary noise, can be regarded as composed of pieces of two deterministic 
equations having two velocity states. It is clear that following this point of 
view, we can perform some extensions: (i) on the number of states of the 
system, (ii) on the laws governing the velocity of the dynamical variable X{t), 
(iii) on the statistics of the generating randomness. Based on these ideas we 
can formalize a piecewise deterministic process as follows: 

(a) The dynamics is described by the equation: 

f^^,(X) (3) 

where As is a function chosen randomly in a set of {Ai, . . . , As} functions. 
Given As, we say that the dynamics is in the state s. We require that 
As{x) be Lipschitz continuous, so that for fixed s, X{t) exists and is 
unique. 

(b) The initial condition is settled either by the Cauchy problem to Eq. (3), 
i.e. X{to) = Xq, and by the initial state of the process. 

(c) States are characterised by a Poisson statistics of switching times: 

l^s e-^^* (4) 

with fis e {/xi, . . .,iis}- 

(d) At every switching time the dynamics can change from the state j (Aj) to 
the state i (Aj) according to a transition probability matrix (or stochastic 
matrix) : 

Qij (5) 
having the following fundamental [2] properties: < q^j < 1 and J2i=i Qij — 
1, for each i,j&l,...,S. 

We can easily visualize this process as a sequence of solutions of Eq. (3), 
where the end point of each is just the initial condition for the successive 
integral curve, so that this process is continuous, because at switching times 
only the first derivative of X{t) is involved. Assumption (c) makes the process 
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be Markovian [1,9], because the exponential statistics, regarded as generated 
by an infinitesimal time process, corresponds to a constant switching rate. 



2. 1 The Liouville- Master Equation 

The statistical description of the process (a,b,c,d) is performed by the Liouville- 
Master Equation for the distribution functions Fs{x,t) — J^^Ps{x' ,t)dx' , for 
each state s = {1, . . . , 5"} of the system: 

s 

dtF,{x, t) + A,{x) d,F,{x, t)^Yl QsjFjix, t), (6) 

where Qsj = {qsj —Ssj)fJ's, and 6sj is the Kronecker's symbol. Eq. (6) is written 
for Fs{x,t) rather than the densities Ps{x,t), but they are connected by a 
derivative. We remind that Ps{x,t) > and Fs{x,t) is non-decreasing in x. 
The Eq. (6) is solved for the Cauchy initial conditions: 

F,{x,0) = Fos{x) = r pos{x')dx' (7) 

J —oo 

given for each state s of the system. Our aim is to find solutions t), ■ ■ • , 
Fs{x, t) of Eq. (6), and then the total distribution function: J^{x, t) = Zlf=i Fs{x, t), 
representing the distribution of the process at time t, regardless of the state 
of the equation. Note that the conservation of probability measure requires 
that both conditions lim,j.^^T{x^t) = 1 and limx^-oo^{x,t) = 0, must be 
satisfied. This means that Eq. (6) is also subject to a limit problem. We note 
that if Ai{x) = —X + 1, A2{x) = —x — 1, /ii = ^2 = /U, Qij = I for i ^ j 
and Qij — ioT i — j, then the filtered random telegraph process (1) and its 
associated Liouville-Master Equation (2) are recovered. 

This system of equation is linear hyperbolic, with non constants coefficients. 
In what follows we focus our attention on solutions J-'{x, t) having an equilib- 
rium, but the numerical scheme can be extended to diffusion processes too. 
Conditions for the existence of equilibrium solution can be conjectured by us- 
ing simple dynamical considerations [17]. If all dynamical equations (3) own 
only attraction points and all of these are contained into the intersection of 
the basin of attraction of each, then a process starting from this region will 
never escape. Whence, there should exists a region fl where the process is 
confined and a stationary distribution J-'eq{x) — \imt^^J^{x,t) exists. 
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3 The numerical solution 

Due to his hyperbolic nature, we use the upwind scheme [3,4] to numerically 
solve Eq. (6). We define the uniform discrete mesh points {xkiU) with steps 
Aa; and At on the domain (i7x [0, T]), so that the coefficient Ai{x) is written as 
^Ak = Ai{xk) and the discrete solution as = Fi{xk,ti). The upwind scheme 
takes the discrete derivative x-direction according to the sign of As{x), so that 
we obtain the following difference discrete scheme for 'F^: 



1/ — 1 corresponds to take in Eq. (8) the spatial discrete right-derivative; u — 
the left one. This is an exphcit scheme of first order, where the solutions 'F^, 
at time step n, are calculated from the Cauchy initial conditions 'F^ = Foi{xk)- 

3.1 Convergence analysis 

Stability conditions for the upwind method are known for linear hyperbolic 
systems as Courant-Friedrichs-Lewy (CFL). If Q = the numerical scheme (8) 
is stable when | Ag At/ Ax |< 1 is fulfilled. If Q 7^ we can use the approach 
of Ref. [3] to state that the discrete operator of (8) is hmited by an exponential 
grow in time, then from the Lax equivalence theorem, the consistency ensures 
convergence. In what follows, convergence of 'F^' to Fi{xk,tn) is proved from 
the global error. It results limited by a linear, rather than exponential, grow 
in time. This remarkable property is a direct consequence of the stochastic 
matrix qij. 

Theorem 1 Let Fi{x,t) E C^'^iQ x [0,T]) be a solution of Eq. (6) with 
\Ai{x)\ < M and Cauchy conditions (7) for x E VL and I G {1, . . . , 5*}, where 
is a closed set of the real axis. Fixed a uniform mesh on (Q x [0,T]), with 
spatial discretization of step-size Ax — Xk+i — x^, and temporal discretization 
At — tn+i — tji'. 



then the global error \\ fe5!||oo,i = S/^i maxfe | fe^i'l = ||F/(xfc,T)— 'F^'Hoo,!, of the 
numerical solution 'F^ resulting from the scheme (8), computed at the time 
step n = T/ At is bounded by \\ 'e^||oo,i < || fe^||oo,i + T£ , where £ = 0{Ax), 
so that 'FjJ* converges to Fi{xk,tn) as Ax — > for {x,t) E fl x [0,T]. 




1 if % < 
if % > 



(8) 



where u — < 
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PROOF. We search for convergence of the numerical method by studying 
Eq. (8) for the global error fe^ — Fi{xk,ti) — '■F^: 



Li+l 



Ax 



+ 1:Qis 'el 



'SI At 



(10) 



where ''£1 is the local truncation error that depends on the second spatial and 
temporal derivative of the solution: 

where and rji are unknown points, and a = At/ Ax. Let ''E'^ = max^ 'e^ 
and = maxfc ''SI , taking modulus and using the triangle inequality gives: 

H+'l < (|l - a\ 'Ak\ + AtQu\ + a\ %{) 'E' + E.^; Qis 'E'+ ^^^^ 

'S' At 

By using the condition of Eq. (9), the first modulus in the right side is positive: 

l-a\%\ + AtQii>0 yi,k (12) 



then we have: 



<{1 + At Qu) 'E' + At E Qis "E' + At 



(13) 



This inequality is valid for all mesh points Xfe and the r.h.s does not depends 
on A;, then we can write: 



lei 



At 



(14) 



Taking the sum over all the states /, we obtain the inequality written for the 
norm on the states: 



E'+'\\i < \\{I + AtQ)\\i\\E% + \\S%At 



(15) 



The matrix I + AtQ gives us the information about the convergence of the 
numerical method (8). By noting that (/ + AtQ)ii > because of (12) and 
Qis = qislJ's >Oiorl^s, we have: 



II (/ + At g)||i = max, El \Sis + At{qis - 6is)fis)\ 

= maxs {J2i Sis{l- AtHs) + At^s J2s Qis) = 1 



(16) 



because the fundamental property (5) of the stochastic transition matrix qij. 
The error of the numerical method at the time step n is limited 



|^"||i < \\E°\\i + nAtS 
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provided that the maximum local truncation error £ — maxj is 0{^x) 
for any {x,t) e O x [0,r]. 



3.2 N on- decreasing property 



The solutions of Eq. (6) are non-decreasing function in space. We search for 
the condition that the numerical solutions 'F^ own the same property, or 
equivalently, that the density probability distributions, defined as h>k — { '^fe ~ 
^F^_i)/Ax will never becomes negative for all i > 0. 

Theorem 2 Let 'F^ the numerical solution obtained via the numerical scheme 
(8), and - ^F^_^ > 0- If the CFL condition (9) is satisfied, then ^F^ - 
H-i > Vi > 0, "iiyk. 



PROOF. By subtracting the solution at meshes point k and A; + 1, we find: 

- %t\ = n - n-i + (-^ 'H+u + ^ H+u-i ^^^^ 

+ ^k+i^-2 + l^s Qls [ ^k- ^k-1 ) ) 



Let Ax 'p^^ = - ^Fltl, we have: 



'Pl^' - ¥k + At + E [Qis - %^^^) k) (18) 

where Qis = {qis - (^zs)a*s, that is: 



(19) 



= (1 - At - Qu)) Vk + At E^^; {Qis %) + 

At I //I I L,j 

^^.1 ^k+u-l\ Pk+2v-l 

We have to ensure that 'f^^ will never becomes negative. Suppose that exists 
i such that 5pa; > for all /c, /, then 

> (1 - At (I /Ax - Q«)) H 

because of Qis > for / 7^ s. Using the CFL restriction for At, we find 
^fe^^ > H ^ 0- Now, being ^° > 0, by induction we find the thesis. 



We observe the CFL condition is a sufficient condition for both convergence 
and non-decreasing of the numerical solution 'F^ (or non negativity of ^\). 
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3.3 Smearing 



In the theorem 1 we have restricted our analysis to solutions having both 
second spatial and temporal derivative, but in general local discontinuities can 
be considered because of the upwind method tends to regularize the solution. 
This is explained in terms of modified equation [4] modeling the numerical 
method. By substituting a function F(a;, t) into Eq. (8), expanding in Taylor's 
series and preserving the second order terms, in matrix notation we get: 

d^F + A{x)d,F + ^dlF - ^^^dlF -QFc^O. (20) 

Taking the derivatives by and df of this equation, we can write: 

d^F ~ A(d,A - Q)d,F + QdtF + A^dlF 
that inserted into Eq. (20) give us: 



(/ - f Q) dtF + A{x) [l + f{d,A - Q)] d^F -QF^ 

A{x)Ax /j _ A{x)At \ Q2p 
V Aa; / a; 



(21) 
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When At — s> we recover Eq. (6) showing the consistency of this modified 
equation. Coefficients of dt and d^x are diagonal matrices and Eq. (21) owns, 
loosely speaking, the diffusive properties of a parabolic equations. The upwind 
method (8) is of the second order with respect to Eq. (21), so that the numer- 
ical solution ^Fl is closer to this modified equation than Eq. (6), showing a 
diffusive character. 



4 Numerical tests 



We perform numerical tests by assigning some values to the parameters model. 
Then we solve the Cauchy problem for the cvolutory equation (6), using the 
upwind method (8), and plot the total density probability p\ = Yl'i=i H- 

In order to test the validity of the numerical scheme versus the CFL condition, 
we study a system of a fluctuation-dissipation process having four relaxation 
states: As{x) = —'jsX+Ws with s = {1, 2, 3, 4}, where the switching times hap- 
pen according to the Poisson processes of Eq. (4) . This model is solved for the 
initial condition (Riemann problem) Fog{x) — H{x)/A (or Posi^) = <^(^)/4), 
up to time T, where H{x) is the Heaviside function. As discussed in Sect. 2 
this kind of system the density distribution always reaches an equilibrium, i.e. 
\imt^^p(x,t) = Peq{x) with X ^ Vt = j^min^ {"y^} ; max^ {^}J- The discon- 
tinuous initial condition does not fulfill the hypothesis for convergence, but. 
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Fig. 1. On the left: convergent solution p{xj^,T) of Eq. (6) when the CFL condition 
(9) is satisfied. On the right: non-convergent solution when that is not satisfied. 

as above mentioned, it is subject to smearing that tends to reguralize disconti- 
nuities (see the third picture of Fig. (2)). In Fig. (1) we show the total density 
distribution function ^(a;^, T) when the temporal step At meet or not the CFL 
condition. We use the following parameters: jis = 4, qis — 0.25, 7s = 10~^, 
Wi = 1, = -1, W3 = 2, = -2. By using Ax = 4.004, we get from the 
CFL condition of Eq. (9) that At < Atmax = 0.250063. On the left side of 
Fig. (1), we used At = 0.225056 < At^ax and T = 545.5. On the right side: 
At = 0.5 > Atmax and the solution is not convergent. In Fig. (2) we plot the 
total probabihty p{xk, U) for fixed times to see how the density distribution 
evolves. For these plots /ij — 0.2, 7^ = 0.1, the remaining parameters are same 
to previous test. These results are compared to histograms obtained by Monte 
Carlo simulation of the dynamical equations (3), (4) and (5). Histograms are 
generated by collecting the final values X{t) of the sample paths. 



5 Summary and conclusions 

In this communication we have introduced the continuous piecewise determin- 
istic Markov processes, as an extention of random telegraph processes, and 
the related Liouville-Master Equation for the distribution functions. The up- 
wind method has been successfully applied for solving this equation. We have 
found a Courant-Priederichs-Lewy condition ensuring both convergence and 
non- decreasing property of the numerical distribution function. The numeri- 
cal solution fits with that we obtained by performing Monte Carlo simulation 
of the process. These results are an introductory work for this general class 
of processes. We trust that the piecewise deterministic processes have a po- 
tential, as models for non-deterministic systems, large range of application 
as confirmed by the recent literature involved on the subject, thus justify- 
ing the development of numerical methods for Eq. (6). Many aspects should 
be further investigated, such as: the asymptotic stability for stiff parameters, 
the behaviour for discontinuous distributions, the development of higher order 
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Fig. 2. Temporal evolution of the density probability distribution functions p{x,t) 
at fixed time. S = 4, At = 0.009, Ax = 0.04004, i^j = 0.2, = 0.1. Prom top left 
to bottom right: t = {0, 4, 8, 12, 16, 20}. Monte Carlo's histograms in dashed lines. 

numerical methods, the extention of the model to non-Markovian and many 
dimensional processes. 
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